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The stability of different phases of the three-dimensional non-relativistic electron gas is analyzed 
using stochastic methods. With decreasing density, we observe a continuous transition from the 
paramagnetic to the ferromagnetic fluid, with an intermediate stability range (25 ± 5 < r s < 35 ± 5) 
for the partially spin-polarized liquid. The freezing transition into a ferromagnetic Wigner crystal 
occurs at r s — 65 ± 10. We discuss the relative stability of different magnetic structures in the solid 
phase, as well as the possibility of disordered phases. 
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Ever since the pioneering work of Wigner and Seitz |l]] 
on the cohesive energy of metals, the calculation of the 
ground state energy of the interacting electron gas (quan- 
tum one component plasma, Fermi gas or simply jellium) 
became the object of considerable theoretical interest [|| . 
Indeed, the electron gas provides the simplest model in 
which non-trivial magnetic structures and electron local- 
ization can be realized by varying a single parameter, 
namely the average electron density p. 

In the present paper we investigate the relative sta- 
bility of various broken symmetry phases of the non- 
relativistic three-dimensional electron gas, both fluid and 
solid, using stochastic methods. We report on a macro- 
scopic instability of the Fermi gas which involves par- 
tial spin-polarization (weak ferromagnetism) , in such a 
way that the paramagnetic to ferromagnetic (full spin- 
polarization) transition is not first order, as previously 
assumed, but a continuous one. Moreover we find that 
the transition to a Wigner crystal occurs at a significantly 
larger density than the value commonly accepted || , and 
that near the quantum freezing transition the fee and bec 
crystal phases are nearly degenerate. 

The jellium model consists of N electrons enclosed in 
a box of volume f2 (periodically repeated in space) in the 
presence of a neutralizing background of positive charge. 
Two parameters characterize its zero temperature phase 
diagram, namely the particle density p — N/Cl and the 
spin polarization £ = | JV-f — iVjJ/iV, where iVjn) is the 
number of spin-up(down) electrons (N = JV-j- +iVj_). The 
system is governed by the Hamiltonian (Hartree a.u.) 
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where and pj) are the position and linear momentum 
of particle i, and A is a constant representing the ef- 
fect of the background. Since we are interested in the 
macroscopic properties of this model system, the ther- 
modynamic limit (N, Q — > oo, keeping p constant) is to 
be performed in the end by finite size extrapolation. 



Elementary scaling arguments indicate that the kinetic 
energy term in (|l|) goes as l/r^ (r s is the Wigner-Seitz 
radius in units of the Bohr radius and its relation to 
the density is p~ x = -j-?~s), while the potential energy 
scales as l/r s . Depending on the relative strength be- 
tween Coulomb and kinetic energies, we can characterize 
three different regimes: the weak (r s < 1), intermediate 
(1 < t s < 10), and strong (r s > 10) Coulomb coupling 
regimes. The random-phase approximation Q] provides 
an accurate description of the weak-coupling regime. The 
intermediate coupling region, of direct interest for den- 
sity functional calculations, has been extensively studied 
by numerical |3|Q] and semi-analytic methods ||. Not 
surprisingly, the least known regime is the strongly cor- 
related one, for which an early quantum Monte Carlo 
calculation || is still the most authoritative study. 

To delve into the strong coupling regime we employ the 
variational (VMC) and diffusion (DMC) quantum Monte 
Carlo methods. The starting point is provided by a vari- 
ational wave function of the Jastrow type: 

$t(TI) = J[R>, S] det T [v3] ■ det;[y>] , 9 T (TZ) G H , 

where detj(j)[<p] is a spin up(down) Slater determinant of 
one-electron orbitals ip that are either plane waves (fluid 
phases) or localized functions (crystal phases). In the 
equation above, 1Z = (ri, • • • , rjv) and E = {o\, ■ ■ ■ ,<tn) 
represent the full set of positions and spins of the elec- 
trons in the system. Considering only 2-body correla- 
tions, the Jastrow factor can be written as: J[TZ, S] = 

Yli<j ex P[ v ij(\ r i ~ r j I)]- F° r an infinite system the op- 
timal Vij(r) should decay as 1/r at large distances to 
account for the correct plasmon dispersion relation. For 
large but finite systems, however, that is not the case 
and to reproduce the long distance behavior without re- 
sorting to expensive evaluations of v's by Ewald sums, 
we consider a finite range Vij , vanishing outside a sphere 
of radius Rt tangential to the unit cell: = if 

r > Rt- We have verified that this truncation involves an 
insignificant loss of correlation energy at the variational 
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level 0. Inside the sphere My, which is different for like 
and unlike spin electrons, is given by: 



and 



i(r) = v tj (r)+A l3 exp[ 
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The {ciij} are determined by the electron-electron cusp 
condition [§, while {6y, Ay , Oij } and are 

variational parameters. The {/&} are chosen in such a 
way that and its two first derivatives are continu- 
ous everywhere. The one-electron orbitals used for the 
crystal phases are exponentials of a Pade function: 



<£j(r) = exp 
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where the (positive) constants k\ , lt2 are also variational 
parameters, and the fixed vectors {Rj, j = l,N} are 
distributed on a regular (bcc or fee) lattice. This an- 
alytic form allows tpj to change from Gaussian close to 
each Rj to a decaying exponential far from the local- 
ization centers. The fluid states are eigenstates of zero 
total momentum, while the solid ones are eigenstates of 
finite lattice translations with vanishing total crystal mo- 
mentum. Both are eigenstates of the z-component of the 
total spin S z with eigenvalue HM = h(N^ — N±)/2. All 
the free coefficients are optimized by the variance mini- 
mization technique introduced in Ref . || . The optimized 
wave functions are used to drive the DMC simulation || . 

Several features of the phase diagram discussed below 
depend upon the accurate evaluation of tiny energy dif- 
ferences. To provide a basis to estimate the reliability 
of our results, we report here the relevant aspects of our 
computation, including an estimate for the statistical and 
systematic errors. First of all, because of the high qual- 
ity of the wave functions, and the relatively long runs, 
statistical errors are the least important source of uncer- 
tainties: for all r s , the statistical error is less than 1 % 
of the correlation energy. Moreover, again because of the 
high variational freedom in the wave function, the energy 
gain in going from VMC to DMC is relatively small. As 
expected, the difference SE = E^J 10 - Ef D f c is largest 
at low r s JlO|, because the n-body contributions not in- 
cluded in our Jastrow function become more important 
at high density jfl). Moreover, for each r s we find that 
8E(Q is larger for £ = than for C, = 1, and, in between, 
the C dependence of SE is well represented by the simple 
interpolation: 5E(Q = SE{0) + (6E(l) - SE(0)) ( 2 . 

For the fluid phase the most important source of uncer- 
tainty in our calculations is the finite-size extrapolation 



to the thermodynamic limit. We applied the extrapola- 
tion scheme proposed in Ref. ||: 

Eeo(C) = 

E N (C) - At N (C) + (N/b(C) - l/AvxiC))- 1 , (5) 

where C = ( r sX)> -E-aKC) is the total energy of the fi- 
nite periodic system, b and Ego are fitting parameters, 
Atjv and Avjy describe the size dependence of Hartrcc- 
Fock kinetic and exchange energies (At^ — t^ F — t^f , 
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). The parameter E^ is the extrap- 



olated value of the total energy per electron. Following 
Ref. |J] we assume b(C) = b (r s ) + bi(r s ) ( 4 . The param- 
eters bo(r s ) and bi(r s ) are obtained from the VMC total 
energies of an extended set of systems (N — 1062, 1450, 
and 1930 for ( = 0, N = 531, 725, 965 for ( = 1). This 
form is used to extrapolate all the energies computed 
by VMC and DMC for unpolarized, partially and fully 
polarized fluid systems. For lower densities the extrapo- 
lation is less critical but still important. We verified that 
for the crystal phases the finite size extrapolation error 
is comparable to the statistical one if N > 500. There- 
fore, the results for N — 686 (bcc) and V = 864 (fee) are 
assumed to be equal to the N — > oo limit. 

Extensive DMC calculations are performed for systems 
with V = 1062 (paramagnetic fluid), N = 725 (ferro- 
magnetic fluid), V = 686 (ferro- and antiferro-magnetic 
bcc crystal), V = 864 (ferromagnetic fee phase). For 
r s < 10 our results agree, within the error bar, with 
those reported in Ref. [Q. Representative values for the 
correlation energy extrapolated to the thermodynamic 
limit are reported in Table I, and Fig. [l] shows a phase 
diagram displaying the stability range of the paramag- 
netic fluid, ferromagnetic fluid, and bcc (ferromagnetic) 
Wigner phases. It is apparent that a transition from the 
paramagnetic to a ferromagnetic fluid phase takes place 
around r s ss 25, while from the fluid to the crystal it is 
at r s ss 65. We found that the fee is slightly less stable 
than the bcc structure, and more stable than the fluid 
phase at r s = 70. The total energy difference between 
the two Wigner phases (10 -3 eV/electron) is five times 
smaller than the one between fluid and bcc crystal, and 
is, therefore, at the limit of our resolution. 

The close competition of the fluid and different crys- 
tal forms for r s « 70 suggests that a metastable amor- 
phous phase could be formed in this region of the phase 
diagram, because the large number of available config- 
urations, and the tiny energy differences could make 
the crystallization kinetics exceedingly slow. Preliminary 
calculations, in which the regular lattice of {Rj} vectors 
arc displaced by random amounts (up to 25 % of the 
nearest neighbor distance), show that the total energy 
increases only very slowly with increasing disorder. 

We now turn to the analysis of the spin dependence 
of the total energy near the ferromagnetic and freezing 
transitions. To ease the analysis, we perform a detailed 
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investigation at the VMC level, and then test the validity 
of the resulting picture by a few DMC calculations. 

Regarding the magnetic phase transition, we compute 
the total energy for several sizes and partial spin polariza- 
tions (l2) in the range 0.8 < r s < 30. The extrapolated 
results for r s = 25 are reported in Fig. |2[ It is clear that, 
at this density, the ground state has partial spin polar- 
ization (weak ferromagnet). The full set of data is used 
to fit an interpolation for the total energy (See Ref. Q): 

( £l (r s )- eo (r s ))-(r 1 (r s )C 2 +r 2 (r s )C 4 ) , (6) 

where E^f (£) is the Hartree-Fock energy. In this ex- 
pression eo(i)(r s ) is the correlation energy of the para- 
magnetic(fully polarized) fluid, and together with Ti(r s ) 
and r2(r s ) are fitted to the functional form 



G{r s ) = -2A(l + air a ) 



In 



2A(p ir l /2 + /3 2 r s + fori' 2 + /? 4 r s 2 ) 



(7) 



whose motivation, and connection to the exact high- 
density limit for e ^ are discussed in Ref. [Q. This 
expression for the energy is minimized with respect to £ 
at fixed r s to identify the ground state spin polarization. 
The stable C, as a function of r s is shown in the inset 
of Fig. H It is manifest that the magnetic transition 
is a continuous one, in contrast to what is predicted on 
the basis of simpler interpolations, like the (widely used) 
form proposed by Perdew and Zunger |l4| . 

We emphasize that the interpolation used to determine 
the ground state polarization is based on VMC calcula- 
tions. Nevertheless, as mentioned above, the total en- 
ergy difference between the VMC and DMC results is 
small, and has a smooth and predictable dependence on 
£. Adding SE(() to the fitted total energies does not 
change the transition from continuous to discontinuous, 
although it moves it slightly towards higher values of r s 
(e.g. £ = 0.5 polarization is stable at r s — 26 in VMC, 
while it is at r s = 30 in DMC). We also remind that our 
calculations rely on the fixed-node (FN) approximation, 
and, at present, we cannot exclude that releasing the 
nodal constraints could change the nature of the mag- 
netic transition. However, recent computations reported 
in Ref. [[llj show that the FN error decreases rapidly with 
decreasing density, and is small already at r s = 20. 

Additional evidence supporting the picture of a mag- 
netic instability and a continuous phase transition is pro- 
vided by the analysis of the radial distribution func- 
tion and structure factor, displaying for r s ss 30 long 
range spin-spin correlations even for nominally unpo- 
larized systems (iVf = iVj..) This is illustrated in 
Fig. where the spin-spin radial distribution function 
(gss( r ) — 2[gtt( r ) ~ 5u( r )D i s plotted for two densities 
close to the instability 15 1. At r s = 20, gss displays the 



expected depletion hole for parallel spins, reflecting the 
preference of the system for spin alternation at short dis- 
tances. At ?' s = 25, instead, the spin correlation is small 
and positive at short range, and oscillating at long range, 
pointing to the formation of magnetic domains with par- 
tial spin polarization. This interpretation in terms of 
domains is confirmed by careful analysis of snapshot con- 
figurations. 

Indirect support to our findings comes also from 
recent experiments performed in doped hexaborides 
(Cai_2;La x B6) 16 1. These authors report a weak fer- 
romagnetic phase at low carrier concentration (r s = 28) 
with an ordered moment corresponding to a partial spin 
polarization of about 10 %. Moreover, the observed Curie 
temperature, which is as high as 600 K, is of the same or- 
der of magnitude as the Fermi energy indicating that this 
represents the natural energy scale of the spin system. 

The magnetic phase diagram near the Wigner phase 
transition is somewhat simpler. For 70 < r s < 100 
the exchange energy, although small in absolute terms, 
is still one order of magnitude larger than the correla- 
tion energy, and therefore it is not surprising that the 
ferromagnetic bec structure is more stable than the anti- 
ferromagnetic spin ordering. However, both at the VMC 
and DMC level, the energy difference between these two 
states is very small (~ 1 % of the correlation energy), 



and only the systematic trend E\ 



Jerro 
tot 



< K 



antiferro 
tot 



for 



70 < r s < 100 gives us some confidence in this result. 
As pointed out in the case of the translational disorder, 
because of entropy and kinetics the spin glass phase could 
be the most relevant one for low density electron systems. 

In conclusion, using stochastic methods we have stud- 
ied the magnetic and freezing quantum transitions of the 
fcrmionic one component plasma. We find that, contrary 
to the commonly accepted picture, the paramagnetic to 
the ferromagnetic fluid transition is a continuous one, oc- 
curring over the 25 < r s < 35 density range, which has 
been approached in experiments by doping highly corre- 
lated solids |l(| . The transition to the bec Wigner crystal 
is first order, it occurs at r s = 65 ± 10 and joins two fully 
polarized spin states. Where available, our results for the 
transition densities disagree substantially from those of a 
previous study || . We think that the disagreement is due 
to the extrapolation of the finite size results to the ther- 
modynamics limit. In our case, the importance of this 
extrapolation has been limited by performing computa- 
tions for systems with up to 2000 electrons, in this way 
reducing the corresponding uncertainties. Finally, we re- 
call that our analysis is based on a selected set of broken 
symmetry states. There exists the possibility of other in- 
stabilities such as the inhomogeneous spin-density-wave 
state or superconducting states which we have not con- 
sidered in the present calculation. 
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TABLE I. Minus the DMC correlation energies (in eV) ex- 
trapolated to the thermodynamic limit. For each £, the first 
row refers to the present results, while the second and third 
rows correspond to Ref. [4] and Ref. [3], respectively. The 
statistical error affects the last decimal digit. 



FIG. 1. Total energy difference (eV) times r s of: (•) the 
paramagnetic and the ferromagnetic fluids; (■) the ferromag- 
netic bec crystal and the ferromagnetic fluid. The statistical 
error bar is comparable to the size of the symbols. The ferro- 
magnetic fluid is stable when both symbols are above zero. 

FIG. 2. Total energy (VMC) as a function of £ for a weak 
ferromagnetic state. The statistical error bar is reported for 
the lowest energy point. The inset shows the equilibrium spin 
polarization r s as a function of £. For £ ~ and £ ~ 1 the re- 
sults depend significantly on the interpolation (see text), and 
therefore they have not been reported in the figure. 

FIG. 3. Spin-spin correlation function (VMC) near the 
magnetic instability. 
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